graph_title = "hogo_persons_WLS_trend")
graph_hogo_persons_WLS_trend_covar
estimates_hogo_persons_WLS_trend_covar <- df_estimates 　 #for robustness check
# Chunk 14
# DID estimation
estimation_results <- dynamic_DID_OLS_notrend(dataset = df_analysis,
outcome_var = df_analysis$yoy_persons_receive,
treat_var = df_analysis$unemploy_shock_diff2)
texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_persons_OLS_notrend")
# Event study graph
graph_yoy_hogo_persons_OLS_notrend <- event_study_graph(data = df_estimates,
graph_title = "yoy_hogo_persons_OLS_notrend")
graph_yoy_hogo_persons_OLS_notrend
estimates_yoy_hogo_persons_OLS_notrend <- df_estimates 　 #for robustness check
# Chunk 15
# DID estimation
estimation_results <- dynamic_DID_WLS_notrend(dataset = df_analysis,
outcome_var = df_analysis$yoy_persons_receive,
treat_var = df_analysis$unemploy_shock_diff2)
texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_persons_WLS_notrend")
# Event study graph
graph_yoy_hogo_persons_WLS_notrend <- event_study_graph(data = df_estimates,
graph_title = "yoy_hogo_persons_WLS_notrend")
graph_yoy_hogo_persons_WLS_notrend
estimates_yoy_hogo_persons_WLS_notrend <- df_estimates 　 #for robustness check
# Chunk 16
# DID estimation
estimation_results <- dynamic_DID_OLS_trend(dataset = df_analysis,
outcome_var = df_analysis$yoy_persons_receive,
treat_var = df_analysis$unemploy_shock_diff2)
texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_persons_OLS_trend")
# Event study graph
graph_yoy_hogo_persons_OLS_trend <- event_study_graph(data = df_estimates,
graph_title = "yoy_hogo_persons_OLS_trend")
graph_yoy_hogo_persons_OLS_trend
estimates_yoy_hogo_persons_OLS_trend <- df_estimates 　 #for robustness check
# Chunk 17
# DID estimation
estimation_results <- dynamic_DID_WLS_trend(dataset = df_analysis,
outcome_var = df_analysis$yoy_persons_receive,
treat_var = df_analysis$unemploy_shock_diff2)
texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_persons_WLS_trend")
# Event study graph
graph_yoy_hogo_persons_WLS_trend <- event_study_graph(data = df_estimates,
graph_title = "yoy_hogo_persons_WLS_trend")
ggplotly(graph_yoy_hogo_persons_WLS_trend)
estimates_yoy_hogo_persons_WLS_trend <- df_estimates 　 #for robustness check
results_yoy_hogo_persons_WLS_trend <- estimation_results # for only-post DID table
# Chunk 18
# DID estimation
estimation_results <- dynamic_onlypost_DID_WLS_trend(dataset = df_analysis,
outcome_var = df_analysis$yoy_persons_receive,
treat_var = df_analysis$unemploy_shock_diff2)
texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_persons_WLS_trend")
# Event study graph
graph_yoy_hogo_persons_WLS_trend_onlypost <- event_study_graph(data = df_estimates ,
graph_title = "yoy_hogo_persons_WLS_trend")
ggplotly(graph_yoy_hogo_persons_WLS_trend_onlypost)
estimates_yoy_hogo_persons_WLS_trend_onlypost <- df_estimates #for robustness check
results_yoy_hogo_persons_WLS_trend_onlypost <- estimation_results # for only-post DID table
# Chunk 19
# DID estimation
estimation_results <- dynamic_DID_OLS_notrend_covar8Xcovid_months(dataset = df_analysis,
outcome_var = df_analysis$yoy_persons_receive,
treat_var = df_analysis$unemploy_shock_diff2)
texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_persons_OLS_notrend")
# Event study graph
graph_yoy_hogo_persons_OLS_notrend_covar <- event_study_graph(data = df_estimates,
graph_title = "yoy_hogo_persons_OLS_notrend")
graph_yoy_hogo_persons_OLS_notrend_covar
estimates_yoy_hogo_persons_OLS_notrend_covar <- df_estimates 　 #for robustness check
# Chunk 20
# DID estimation
estimation_results <- dynamic_DID_WLS_notrend_covar8Xcovid_months(dataset = df_analysis,
outcome_var = df_analysis$yoy_persons_receive,
treat_var = df_analysis$unemploy_shock_diff2)
texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_persons_WLS_notrend")
# Event study graph
graph_yoy_hogo_persons_WLS_notrend_covar <- event_study_graph(data = df_estimates,
graph_title = "yoy_hogo_persons_WLS_notrend")
graph_yoy_hogo_persons_WLS_notrend_covar
estimates_yoy_hogo_persons_WLS_notrend_covar <- df_estimates 　 #for robustness check
# Chunk 21
# DID estimation
estimation_results <- dynamic_DID_OLS_trend_covar8Xcovid_months(dataset = df_analysis,
outcome_var = df_analysis$yoy_persons_receive,
treat_var = df_analysis$unemploy_shock_diff2)
texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_persons_OLS_trend")
# Event study graph
graph_yoy_hogo_persons_OLS_trend_covar <- event_study_graph(data = df_estimates,
graph_title = "yoy_hogo_persons_OLS_trend")
graph_yoy_hogo_persons_OLS_trend_covar
estimates_yoy_hogo_persons_OLS_trend_covar <- df_estimates 　 #for robustness check
# DID estimation
estimation_results <- dynamic_DID_WLS_trend_covar8Xcovid_months(dataset = df_analysis,
outcome_var = df_analysis$yoy_persons_receive,
treat_var = df_analysis$unemploy_shock_diff2)
texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_persons_WLS_trend")
# Event study graph
graph_yoy_hogo_persons_WLS_trend_covar <- event_study_graph(data = df_estimates,
graph_title = "yoy_hogo_persons_WLS_trend")
graph_yoy_hogo_persons_WLS_trend_covar
estimates_yoy_hogo_persons_WLS_trend_covar <- df_estimates 　 #for robustness check
results_yoy_hogo_persons_WLS_trend_covar <- estimation_results # for only-post DID table
ggplotly(graph_yoy_hogo_persons_WLS_trend_covar)
# DID estimation
estimation_results <- dynamic_onlypost_DID_WLS_trend_covar8Xcovid_months(dataset = df_analysis,
outcome_var = df_analysis$yoy_persons_receive,
treat_var = df_analysis$unemploy_shock_diff2)
#texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_persons_WLS_trend")
# Event study graph
graph_yoy_hogo_persons_WLS_trend_covar_onlypost <- event_study_graph(data = df_estimates ,
graph_title = "yoy_hogo_persons_WLS_trend")
ggplotly(graph_yoy_hogo_persons_WLS_trend_covar_onlypost)
estimates_yoy_hogo_persons_WLS_trend_covar_onlypost <- df_estimates #for robustness check
results_yoy_hogo_persons_WLS_trend_covar_onlypost <- estimation_results # for only-post DID table
# DID estimation
estimation_results <- dynamic_onlypost_DID_WLS_trend_covar8Xcovid_months(dataset = df_analysis,
outcome_var = df_analysis$yoy_persons_receive,
treat_var = df_analysis$unemploy_shock_diff2)
texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_persons_WLS_trend")
# Event study graph
graph_yoy_hogo_persons_WLS_trend_covar_onlypost <- event_study_graph(data = df_estimates ,
graph_title = "yoy_hogo_persons_WLS_trend")
ggplotly(graph_yoy_hogo_persons_WLS_trend_covar_onlypost)
estimates_yoy_hogo_persons_WLS_trend_covar_onlypost <- df_estimates #for robustness check
results_yoy_hogo_persons_WLS_trend_covar_onlypost <- estimation_results # for only-post DID table
# DID estimation
estimation_results <- dynamic_DID_WLS_trend(dataset = df_analysis,
outcome_var = df_analysis$yoy_households_receive,
treat_var = df_analysis$unemploy_shock_diff2)
texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_households_WLS_trend")
# Event study graph
graph_yoy_hogo_households_WLS_trend <- event_study_graph(data = df_estimates,
graph_title = "yoy_hogo_households_WLS_trend")
graph_yoy_hogo_households_WLS_trend
estimates_yoy_hogo_households_WLS_trend <- df_estimates 　 #for robustness check
results_yoy_hogo_households_WLS_trend <- estimation_results # for only-post DID table
ggplotly(graph_yoy_hogo_households_WLS_trend)
# DID estimation
estimation_results <- dynamic_onlypost_DID_WLS_trend(dataset = df_analysis,
outcome_var = df_analysis$yoy_households_receive,
treat_var = df_analysis$unemploy_shock_diff2)
#texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_households_WLS_trend")
# Event study graph
graph_yoy_hogo_households_WLS_trend_onlypost <- event_study_graph(data = df_estimates ,
graph_title = "yoy_hogo_households_WLS_trend")
ggplotly(graph_yoy_hogo_households_WLS_trend_onlypost)
estimates_yoy_hogo_households_WLS_trend_onlypost <- df_estimates #for robustness check
results_yoy_hogo_households_WLS_trend_onlypost <- estimation_results # for only-post DID table
# DID estimation
estimation_results <- dynamic_onlypost_DID_WLS_trend(dataset = df_analysis,
outcome_var = df_analysis$yoy_households_receive,
treat_var = df_analysis$unemploy_shock_diff2)
texreg::screenreg(l = estimation_results , include.ci = FALSE, digits=3,)
# DID estimates and 90% CI
df_estimates <- DID_coefficients_main(DID_results = estimation_results,
treat_var = "unemploy_shock_diff2",
estimation_label = "yoy_hogo_households_WLS_trend")
# Event study graph
graph_yoy_hogo_households_WLS_trend_onlypost <- event_study_graph(data = df_estimates ,
graph_title = "yoy_hogo_households_WLS_trend")
ggplotly(graph_yoy_hogo_households_WLS_trend_onlypost)
estimates_yoy_hogo_households_WLS_trend_onlypost <- df_estimates #for robustness check
results_yoy_hogo_households_WLS_trend_onlypost <- estimation_results # for only-post DID table
# Chunk 1: setup
Sys.setenv(LANG = "en") #English
knitr::opts_chunk$set(echo = TRUE)
# Chunk 2
rm(list = ls())
path <- getwd()
setwd(path)
# packages
pacman::p_load(tidyverse, patchwork, plotly)
# Font for windows and mac
if (stringr::str_detect(path, pattern="Users")){
theme_set(theme_classic(base_size = 10, base_family = "HiraginoSans-W3"))  # For Mac OS
} else{
theme_set(theme_classic(base_size = 10, base_family = "Arial")) # For Windows
}
options(scipen=10)
# annotation size
annotate_size <- 3
# Chunk 3
# 社会人口統計体系
#統計でみる都道府県のすがた2021
# A 人口・世帯、C　経済基盤、F　労働のデータ
#https://www.e-stat.go.jp/stat-search/files?page=1&layout=datalist&toukei=00200502&tstat=000001149949&cycle=0&year=20210&month=0&tclass1=000001149950&tclass2val=0
#統計でみる都道府県のすがた2021
# 居住
#https://www.e-stat.go.jp/stat-search/files?page=1&layout=datalist&toukei=00200502&tstat=000001149949&cycle=0&tclass1=000001149950&cycle_facet=cycle&tclass2val=0
#社会人口統計体系 (2021)
#Population_per_1_km_2_of_inhabitable_area  = 総人口/可住地面積
## 可住地面積　＝　総面積−森林面積
### 可住地面積の収書は以下のいずれか
### 全国都道府県市区町村調・世界農林業センサス・農林業センサス
#Ratio_of_aged_population 人口推計 2019年10月1日現在人口推計
#Total_population 人口推計 2019年10月1日現在人口推計
# エクセルでデータを編集
Social_Demographic_Stat <- readxl::read_excel("input/social_demographic_stat.xlsx")
# Chunk 4
census_Industrial_Classification_Population <- readr::read_csv("input/census_Industrial_Classification_Population.csv", locale=locale(encoding="CP932"))
# Chunk 5
census_Industrial_Classification_Population <- census_Industrial_Classification_Population %>%
dplyr::rename("prefec_kanji"="地域2015",
"Total_pop"="総数",
"Manufacturing" = "Ｅ製造業",
"Wholesale_Retail_trade"="Ｉ卸売業，小売業",
"Accomodations_eating_drinking_services" = "Ｍ宿泊業，飲食サービス業" ,
"Living_related_personal_services_amusement_services" ="Ｎ生活関連サービス業，娯楽業",
"Secondary_industry"="（再掲）第２次産業" ,
"Tertiary_industry"="（再掲）第３次産業" )
census_Industrial_Classification_Population <- census_Industrial_Classification_Population %>%
dplyr::select(prefec_kanji, Total_pop, Manufacturing, Wholesale_Retail_trade, Accomodations_eating_drinking_services, Living_related_personal_services_amusement_services, Secondary_industry, Tertiary_industry)
# Chunk 6
census_Industrial_Classification_Population <- census_Industrial_Classification_Population %>%
dplyr::mutate(Manufacturing_ratio = Manufacturing/Total_pop) %>%
dplyr::mutate(Wholesale_Retail_trade_ratio = Wholesale_Retail_trade/Total_pop) %>%
dplyr::mutate(Accomodations_eating_drinking_services_ratio = Accomodations_eating_drinking_services/Total_pop) %>%
dplyr::mutate(Living_related_personal_services_amusement_services_ratio = Living_related_personal_services_amusement_services/Total_pop) %>%
dplyr::mutate(Secondary_industry_ratio = Secondary_industry/Total_pop) %>%
dplyr::mutate(Tertiary_industry_ratio = Tertiary_industry/Total_pop)
# Chunk 7
#国勢調査(2015)
#Secondary_industry_ratio ＝（再掲）第２次産業/総数 (単位人)
#Tertiary_industry_ratio＝（再掲）第３次産業/総数 (単位人)
employment_status_survey <- readr::read_csv("input/employment_status_survey.csv",locale=locale(encoding="CP932"))
# Chunk 8
for(i in 1:ncol(employment_status_survey)){
colnames(employment_status_survey)[i] <-  paste(employment_status_survey[[1,i]], colnames(employment_status_survey)[i], sep="_")}
employment_status_survey <- employment_status_survey[-1,]
View(employment_status_survey)
# Chunk 1: setup
Sys.setenv(LANG = "en") #English
knitr::opts_chunk$set(echo = TRUE)
# Chunk 2
rm(list = ls())
path <- getwd()
setwd(path)
# packages
pacman::p_load(tidyverse, patchwork, plotly)
# Font for windows and mac
if (stringr::str_detect(path, pattern="Users")){
theme_set(theme_classic(base_size = 10, base_family = "HiraginoSans-W3"))  # For Mac OS
} else{
theme_set(theme_classic(base_size = 10, base_family = "Arial")) # For Windows
}
options(scipen=10)
# annotation size
annotate_size <- 3
# Chunk 3
# 社会人口統計体系
#統計でみる都道府県のすがた2021
# A 人口・世帯、C　経済基盤、F　労働のデータ
#https://www.e-stat.go.jp/stat-search/files?page=1&layout=datalist&toukei=00200502&tstat=000001149949&cycle=0&year=20210&month=0&tclass1=000001149950&tclass2val=0
#統計でみる都道府県のすがた2021
# 居住
#https://www.e-stat.go.jp/stat-search/files?page=1&layout=datalist&toukei=00200502&tstat=000001149949&cycle=0&tclass1=000001149950&cycle_facet=cycle&tclass2val=0
#社会人口統計体系 (2021)
#Population_per_1_km_2_of_inhabitable_area  = 総人口/可住地面積
## 可住地面積　＝　総面積−森林面積
### 可住地面積の収書は以下のいずれか
### 全国都道府県市区町村調・世界農林業センサス・農林業センサス
#Ratio_of_aged_population 人口推計 2019年10月1日現在人口推計
#Total_population 人口推計 2019年10月1日現在人口推計
# エクセルでデータを編集
Social_Demographic_Stat <- readxl::read_excel("input/social_demographic_stat.xlsx")
# Chunk 4
census_Industrial_Classification_Population <- readr::read_csv("input/census_Industrial_Classification_Population.csv", locale=locale(encoding="CP932"))
# Chunk 5
census_Industrial_Classification_Population <- census_Industrial_Classification_Population %>%
dplyr::rename("prefec_kanji"="地域2015",
"Total_pop"="総数",
"Manufacturing" = "Ｅ製造業",
"Wholesale_Retail_trade"="Ｉ卸売業，小売業",
"Accomodations_eating_drinking_services" = "Ｍ宿泊業，飲食サービス業" ,
"Living_related_personal_services_amusement_services" ="Ｎ生活関連サービス業，娯楽業",
"Secondary_industry"="（再掲）第２次産業" ,
"Tertiary_industry"="（再掲）第３次産業" )
census_Industrial_Classification_Population <- census_Industrial_Classification_Population %>%
dplyr::select(prefec_kanji, Total_pop, Manufacturing, Wholesale_Retail_trade, Accomodations_eating_drinking_services, Living_related_personal_services_amusement_services, Secondary_industry, Tertiary_industry)
# Chunk 6
census_Industrial_Classification_Population <- census_Industrial_Classification_Population %>%
dplyr::mutate(Manufacturing_ratio = Manufacturing/Total_pop) %>%
dplyr::mutate(Wholesale_Retail_trade_ratio = Wholesale_Retail_trade/Total_pop) %>%
dplyr::mutate(Accomodations_eating_drinking_services_ratio = Accomodations_eating_drinking_services/Total_pop) %>%
dplyr::mutate(Living_related_personal_services_amusement_services_ratio = Living_related_personal_services_amusement_services/Total_pop) %>%
dplyr::mutate(Secondary_industry_ratio = Secondary_industry/Total_pop) %>%
dplyr::mutate(Tertiary_industry_ratio = Tertiary_industry/Total_pop)
# Chunk 7
#国勢調査(2015)
#Secondary_industry_ratio ＝（再掲）第２次産業/総数 (単位人)
#Tertiary_industry_ratio＝（再掲）第３次産業/総数 (単位人)
employment_status_survey <- readr::read_csv("input/employment_status_survey.csv",locale=locale(encoding="CP932"))
# Chunk 8
for(i in 1:ncol(employment_status_survey)){
colnames(employment_status_survey)[i] <-  paste(employment_status_survey[[1,i]], colnames(employment_status_survey)[i], sep="_")}
employment_status_survey <- employment_status_survey[-1,]
# Chunk 1: setup
Sys.setenv(LANG = "en") #English
knitr::opts_chunk$set(echo = TRUE)
# Chunk 2
rm(list = ls())
path <- getwd()
setwd(path)
# packages
pacman::p_load(tidyverse, patchwork, plotly)
# Font for windows and mac
if (stringr::str_detect(path, pattern="Users")){
theme_set(theme_classic(base_size = 10, base_family = "HiraginoSans-W3"))  # For Mac OS
} else{
theme_set(theme_classic(base_size = 10, base_family = "Arial")) # For Windows
}
options(scipen=10)
# annotation size
annotate_size <- 3
# Chunk 3
# 社会人口統計体系
#統計でみる都道府県のすがた2021
# A 人口・世帯、C　経済基盤、F　労働のデータ
#https://www.e-stat.go.jp/stat-search/files?page=1&layout=datalist&toukei=00200502&tstat=000001149949&cycle=0&year=20210&month=0&tclass1=000001149950&tclass2val=0
#統計でみる都道府県のすがた2021
# 居住
#https://www.e-stat.go.jp/stat-search/files?page=1&layout=datalist&toukei=00200502&tstat=000001149949&cycle=0&tclass1=000001149950&cycle_facet=cycle&tclass2val=0
#社会人口統計体系 (2021)
#Population_per_1_km_2_of_inhabitable_area  = 総人口/可住地面積
## 可住地面積　＝　総面積−森林面積
### 可住地面積の収書は以下のいずれか
### 全国都道府県市区町村調・世界農林業センサス・農林業センサス
#Ratio_of_aged_population 人口推計 2019年10月1日現在人口推計
#Total_population 人口推計 2019年10月1日現在人口推計
# エクセルでデータを編集
Social_Demographic_Stat <- readxl::read_excel("input/social_demographic_stat.xlsx")
# Chunk 4
census_Industrial_Classification_Population <- readr::read_csv("input/census_Industrial_Classification_Population.csv", locale=locale(encoding="CP932"))
# Chunk 5
census_Industrial_Classification_Population <- census_Industrial_Classification_Population %>%
dplyr::rename("prefec_kanji"="地域2015",
"Total_pop"="総数",
"Manufacturing" = "Ｅ製造業",
"Wholesale_Retail_trade"="Ｉ卸売業，小売業",
"Accomodations_eating_drinking_services" = "Ｍ宿泊業，飲食サービス業" ,
"Living_related_personal_services_amusement_services" ="Ｎ生活関連サービス業，娯楽業",
"Secondary_industry"="（再掲）第２次産業" ,
"Tertiary_industry"="（再掲）第３次産業" )
census_Industrial_Classification_Population <- census_Industrial_Classification_Population %>%
dplyr::select(prefec_kanji, Total_pop, Manufacturing, Wholesale_Retail_trade, Accomodations_eating_drinking_services, Living_related_personal_services_amusement_services, Secondary_industry, Tertiary_industry)
# Chunk 6
census_Industrial_Classification_Population <- census_Industrial_Classification_Population %>%
dplyr::mutate(Manufacturing_ratio = Manufacturing/Total_pop) %>%
dplyr::mutate(Wholesale_Retail_trade_ratio = Wholesale_Retail_trade/Total_pop) %>%
dplyr::mutate(Accomodations_eating_drinking_services_ratio = Accomodations_eating_drinking_services/Total_pop) %>%
dplyr::mutate(Living_related_personal_services_amusement_services_ratio = Living_related_personal_services_amusement_services/Total_pop) %>%
dplyr::mutate(Secondary_industry_ratio = Secondary_industry/Total_pop) %>%
dplyr::mutate(Tertiary_industry_ratio = Tertiary_industry/Total_pop)
# Chunk 7
#国勢調査(2015)
#Secondary_industry_ratio ＝（再掲）第２次産業/総数 (単位人)
#Tertiary_industry_ratio＝（再掲）第３次産業/総数 (単位人)
employment_status_survey <- readr::read_csv("input/employment_status_survey.csv",locale=locale(encoding="CP932"))
# Chunk 8
for(i in 1:ncol(employment_status_survey)){
colnames(employment_status_survey)[i] <-  paste(employment_status_survey[[1,i]], colnames(employment_status_survey)[i], sep="_")}
employment_status_survey <- employment_status_survey[-1,]
View(employment_status_survey)
colnames(employment_status_survey)
employment_status_survey <- readr::read_csv("input/employment_status_survey.csv",locale=locale(encoding="CP932"))
employment_status_survey <- readr::read_csv("input/employment_status_survey.csv",locale=locale(encoding="CP932"))
employment_status_survey <- readr::read_csv("input/employment_status_survey.csv",locale=locale(encoding="CP932"))
employment_status_survey <- readr::read_csv("input/employment_status_survey.csv",locale=locale(encoding="CP932"))
rm(list = ls())
path <- getwd()
setwd(path)
# packages
pacman::p_load(tidyverse, plotly,readxl,DT, lubridate, extrafont)
# Font for windows and mac
if (stringr::str_detect(path, pattern="Users")){
theme_set(theme_classic(base_size = 10, base_family = "HiraginoSans-W3"))  # For Mac OS
base_family <- "HiraginoSans-W3"
} else{
theme_set(theme_classic(base_size = 10, base_family = "Arial")) # For Windows
base_family <- "Arial"
}
newly_confirmed_cases_daily <- readr::read_csv("input/newly_confirmed_cases_daily.csv",col_types = cols(Date= col_date(format="%Y/%m/%d")))
confirmed_cases_cumulative_daily <-readr::read_csv("input/confirmed_cases_cumulative_daily.csv",col_types = cols(Date= col_date(format="%Y/%m/%d")))
deaths_cumulative_daily <- readr::read_csv("input/deaths_cumulative_daily.csv",col_types = cols(Date= col_date(format="%Y/%m/%d")))
JP_Mobility_Report <- readr::read_csv("input/2020_JP_Region_Mobility_Report.csv")
prefec_id <- read_excel("input/prefec_id.xlsx")
prefec_id <- read_excel("input/prefec_id.xlsx")
# Data on the number of cumulative deaths at the end of June 2020
cum_deaths_20200630 <- deaths_cumulative_daily %>%
mutate(year = year(Date),month=month(Date),day=day(Date)) %>%
filter(Date == "2020-06-30") %>%
mutate(.,month = 6,year=2020)
# Recode "prefecture" to "prerec" and "All" to "zenkoku"/後のleft_joinのためにprefectureをprerecに、Allをzenkokuにrecode
cum_deaths_20200630<- cum_deaths_20200630 %>%
mutate(.,prefec = recode(Prefecture, ALL = "zenkoku"))
# cumulative confirmed cases
cum_confirmed_cases_20200630 <- confirmed_cases_cumulative_daily %>%
mutate(year = year(Date),month=month(Date),day=day(Date)) %>%
filter(Date == "2020-06-30") %>%
mutate(.,month = 6,year=2020)
# Recode "prefecture" to "prerec" and "All" to "zenkoku"/後のleft_joinのためにprefectureをprerecに、Allをzenkokuにrecode
cum_confirmed_cases_20200630<- cum_confirmed_cases_20200630 %>%
mutate(.,prefec = recode(Prefecture, ALL = "zenkoku"))
# Data on the number of cumulative deaths at the end of may 2020
cum_deaths_20200531 <- deaths_cumulative_daily %>%
mutate(year = year(Date),month=month(Date),day=day(Date)) %>%
filter(Date == "2020-05-31") %>%
mutate(.,month = 5,year=2020)
# Recode "prefecture" to "prerec" and "All" to "zenkoku"/後のleft_joinのためにprefectureをprerecに、Allをzenkokuにrecode
cum_deaths_20200531<- cum_deaths_20200531 %>%
mutate(.,prefec = recode(Prefecture, ALL = "zenkoku"))
# cumulative confirmed cases
cum_confirmed_cases_20200531 <- confirmed_cases_cumulative_daily %>%
mutate(year = year(Date),month=month(Date),day=day(Date)) %>%
filter(Date == "2020-05-31") %>%
mutate(.,month = 5,year=2020)
#  Recode "prefecture" to "prerec" and "All" to "zenkoku後のleft_joinのためにprefectureをprerecに、Allをzenkokuにrecode
cum_confirmed_cases_20200531<- cum_confirmed_cases_20200531 %>%
mutate(.,prefec = recode(Prefecture, ALL = "zenkoku"))
google_mobility_Japan2020  <- JP_Mobility_Report %>%
mutate(year = year(date),month=month(date),day=day(date)) %>%
# 2020年にデータを限定
filter(year==2020) %>%
# 日本全国のsub_region1がnaなのでreplace_na
replace_na(list(sub_region_1="zenkoku")) %>%
# sub_region1と月次で集計
group_by(sub_region_1,month) %>%
summarise(
retail_and_recreation = mean(retail_and_recreation_percent_change_from_baseline),
grocery_and_pharmacy = mean(grocery_and_pharmacy_percent_change_from_baseline),
parks = mean(parks_percent_change_from_baseline),
transit_stations = mean(transit_stations_percent_change_from_baseline),
workplaces = mean(workplaces_percent_change_from_baseline),
residential = mean(residential_percent_change_from_baseline)
)
# make mobility index
#  year=2020代入、sub_region1をprefecに変換。
google_mobility_Japan2020 <-google_mobility_Japan2020 %>%
mutate(.,prefec=sub_region_1,
year=2020,
google_mobility_index_4vari_average = (retail_and_recreation+grocery_and_pharmacy+transit_stations+workplaces)/4,
mobility_index_6vari_average =(retail_and_recreation+grocery_and_pharmacy+transit_stations+workplaces+parks+residential)/6)
covariates2020may <-left_join(cum_confirmed_cases_20200531, cum_deaths_20200531,by = c("prefec","year","month","day","Date","Prefecture"))
covariates2020may <- left_join(google_mobility_Japan2020, covariates2020may,by = c("prefec","month","year"))
covariates2020may <- left_join(covariates2020may,prefec_id,by="prefec")
covariates <-left_join(cum_confirmed_cases_20200630, cum_deaths_20200630,by = c("prefec","year","month","day","Date","Prefecture"))
covariates2020 <- left_join(google_mobility_Japan2020, covariates,by = c("prefec","month","year"))
covariates2020 <- left_join(covariates2020,prefec_id,by="prefec")
covariates202005 <- covariates2020may %>%
# sub_region_1をドロップするために、ungroupを使う。googlemobliity indexの箇所でgroup_byを使っているため。
filter(month==5) %>%
ungroup() %>%
select(-c(sub_region_1,Prefecture ))
covariates202005 <- covariates202005 %>%
rename(confirmed_cases_cumulative="Confirmed cases(Cumulative)",deaths_cumulative="Deaths(Cumulative)")
renv::init()
renv::snapshot()
renv::restore()
